Statistical Mechanics of Time Independent Non-Dissipative Nonequilibrium States 



Stephen R. William^] and Denis J. Evanfj] 

Research School of Chemistry, Australian National University, Canberra, ACT 0200, Australia 

(Dated: 1st February 2008) 

We examine the question of whether the formal expressions of equilibrium statistical mechanics 
can be applied to time independent non-dissipative systems that are not in true thermodynamic 
equilibrium and are nonergodic. By assuming the phase space may be divided into time independent, 
locally ergodic domains, we argue that within such domains the relative probabilities of microstates 
are given by the standard Boltzmann weights. In contrast to previous energy landscape treatments, 
that have been developed specifically for the glass transition, we do not impose an a priori knowledge 
of the inter-domain population distribution. Assuming that these domains are robust with respect to 
small changes in thermodynamic state variables we derive a variety of fluctuation formulae for these 
systems. We verify our theoretical results using molecular dynamics simulations on a model glass 
forming system. Non-equilibrium Transient Fluctuation Relations are derived for the fluctuations 
resulting from a sudden finite change to the system's temperature or pressure and these are shown to 
be consistent with the simulation results. The necessary and sufficient conditions for these relations 
to be valid are that the domains are internally populated by Boltzmann statistics and that the 
domains are robust. The Transient Fluctuation Relations thus provide an independent quantitative 
justification for the assumptions used in our statistical mechanical treatment of these systems. 

I. INTRODUCTION 

The formal expressions of equilibrium statistical mechanics strictly apply only to ergodic systems that are in 
thermodynamic equilibrium. Thus these expressions only strictly apply to systems which are at the global free energy 
minimum given the system Hamiltonian and the macroscopic thermodynamic state variables (number of particles, 
temperature and pressure or density). For such systems Gibbsian equilibrium statistical mechanics provides an exact 
prescription for how to calculate the various thermodynamic quantities^. However, these prescriptions are routinely 
applied to systems that are not in true thermodynamic equilibrium (for example to metastable liquids 2 , glasses 3 , 
polymorphs 4 and allotropes) . It is often observed empirically that within experimental uncertainties many expressions 
for thermodynamic quantities yield consistent results. In the present paper we provide arguments for why many of 
the results of equilibrium statistical mechanics can be applied to such time independent nondissipative nonequilibrium 
systems. We also point out some of the limits inherent in the application of the formulae of equilibrium statistical 
mechanics to such systems. 

We choose to study the isothermal isobaric ensemble 5 - (externally regulated pressure and temperature) . The methods 
and reasoning we use here can be directly transferred to other ensembles such as the canonical (fixed volume and 
externally regulated temperature) . The Gibbs free energy G, which is the thermodynamic potential for the isothermal 
isobaric ensemble, is related to the partition function A by the equation 

G{N,P ,T) = -k B ThxA(N,P ,T), (1) 

and the partition function is given by the integral 

A = J J dV dTexp[-P{H {T) + P V)}, (2) 

where T = (q,p) is the phase space vector describing the coordinates q and momenta p, of all the N particles in 
the system, Pq is the thermodynamic pressure, and = 1/fcgT where ks is Boltzmann's constant and T is the 
temperature. The integration domain D provides limits for both integrals and extends over all the available phase 
space (r, V). This is ±oo for every component of the generalized momentum, — > oo for the volume V, and over the 
volume for the Cartesian coordinates of the particles. Since the system Hamiltonian Hq(T, V) is single valued, so too 
is the partition function and in turn the free energy. 

If we require the distribution function of a single thermodynamic phase it is necessary that other phases do not 
contribute significantly to the partition function. The full integration domain D may include states that are charac- 
teristic of crystalline states or fluids states. In the thermodynamic limit this does not cause problems because, as we 
shall see, the partition function will be completely dominated by those microscopic domains that have the lowest free 
energy. However the application of these formulae to allotropes or metastable systems does present a problem. The 
standard equilibrium statistical mechanical expressions for variables such as the enthalpy I, the average volume (V) 
and second order quantities such as the specific heat at constant pressure cp may all be computed from a knowledge 
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of the partition function Eq. [2] or equivalently the thermodynamic potential Eq. [TJ If other phases of lower free 
energy exist this computation (from Eq. [2] as written) will strictly speaking be incorrect. 

It is well known that the formulae for thermal properties such as entropy, free energy, temperature and specific heat 
do not hold for dissipative nonequilibrium systems outside the linear response regime^ 7 -. In this paper we examine 
the question of whether they are correct for any nondissipative nonequilibrium systems such as allotropes, metastable 
systems or history dependent glasses. We provide a statistical mechanical theory of time independent, nondissipative, 
nonequilibrium systems. The theory is based on the fact that these systems are nonergodic and individual sample 
systems comprise ergodic domains that do not span all of phase space. We show that if these domains are robust 
with respect to small changes in thermodynamic state variables, a successful statistical mechanical treatment of 
these nonequilibrium systems can be given. We provide direct evidence, from molecular dynamics simulations on a 
model glass former, that the resulting statistical mechanical formulae are satisfied within empirical errors. Finally 
we provide an independent test of the two key elements of our theory: Boltzmann weights within the phase space 
domains and the robustness of those domains. It happens that these two elements are the necessary and sufficient 
conditions for the application of the transient fluctuation relation to finite thermodynamic quenches (in temperature 
or pressure) for such systems^. While the application of thermodynamics to a single time averaged system is usually 
straightforward the application to a ensemble, whose members may be locked in different phase space domains, can 
require modification to the standard formulae. 

In the case of glasses our treatment has some similarities with the energy landscape approach of Stillinger and 
Weber— However there are significant differences; our treatment makes no reference to the inherent structure 
and imposes no a priori knowledge of the inter-domain relative population levels. The energy landscape approach has 
been extended to account for the phenomena of ageing or history dependence by the addition of a Active parameter—. 
Sciortino has convincingly shown that the addition of a single Active parameter is inadequate to deal with glasses, 
which may have different properties at the same temperature and pressure if they are prepared by a different protocol 
(different history dependence)^ and poses the challenge to recover a thermodynamic description "by decomposing the 
ageing system into a collection of substates". The treatment we present here succeeds in doing just that by providing 
a rigorous development of equilibrium statistical mechanics and thermodynamics for ensembles of systems where the 
phase space breaks up into ensembles of domains whose inter-domain dynamics is nonergodic and whose inter-domain 
population levels may not be Boltzmann weighted. 



II. CONDITIONS FOR EQUILIBRIUM 

A dynamical system in equilibrium has the properties that it is nondissipative and that its macroscopic properties 
are time independent. Thus the N-particle phase space distribution function, f(T, V, t), must be a time independent 
solution to the Liouville equation 7 -, 

^/(r', t) = f • v/(r', t) - f(V, t)A = o, (3) 

where A is the phase space compression factor— obtained by taking the divergence of the equations of motion (see 
Eq. [5]) and r' is the extended phase space vector which consists of T and may include additional dynamical variables 
such as the volume V. Since the system is assumed to be nondissipative both the ensemble average (A) and the 
time average A of the phase space compression factor (which is directly proportional to the rate at which heat is 
exchanged with the fictitious thermostat) are zero. The time independent solution to Eq. [3] depends on the details 
of the equations of motion. Equilibrium solutions to Eq. [3] for the equations of motion, suitable for use in molecular 
dynamics simulations, are compatible with Gibbsian equilibrium statistical mechanics 7 -. 

Microscopic expressions for mechanical properties like the pressure, the internal energy, the enthalpy and the volume 
can be derived without reference to Gibbsian statistical mechanics and indeed can be proved to hold for nonequilibrium 
systems including nonequilibrium dissipative systems. 

There are two ways in which the formulae derived from Gibbsian equilibrium statistical mechanics can break 
down. The most obvious way is that the relative weights of microstates may be non-Boltzmann and the exponential 
factor exp[— [3Hq(T)], may be replaced by some other function (either the exponential function itself may be modified 
as in Tsallis statistics^ or the Hamiltonian may be modified to some new function Ho(T) — » B(T,t)Ho(T)). In 
either circumstance the standard expressions for the thermal quantities derived from equilibrium Gibbsian statistical 
mechanics will not be valid. This certainly happens in dissipative nonequilibrium systems where the distribution 
function is not a time independent solution to Eq. 03 

In deterministic nonequilibrium steady states the phase space may break down into ergodically separated domains 
(Each of which will be fractal and of lower dimension than the ostensible phase space dimension. This is a conse- 
quence of dissipation.) However for these steady states, the domains are always exquisitely sensitive to macroscopic 
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thermodynamic parameters since they are strange fractal attractors^. Often a deterministic nonequilibrium steady 
state approaches a unique fractal attractor. As time progresses the distribution function collapses ever closer to (but 
never reaching) the steady state attractor. 

The second way that these expressions may fail is that the system may become nonergodic. In this case three things 
happen, a) Most obviously time averages no longer equal full (domain D) ensemble averages, b) If we take an initial 
microstate the subsequent phase space trajectory will span some phase space domain D a where the initial phase is 
labeled T' a (0). In this case for nondissipative nonequilibrium systems where the domains are robust (i.e. small changes 
in thermodynamic state parameters, to leading order do not change the domain) the standard equations of equilibrium 
statistical mechanics may continue to be valid but in a slightly modified form. We will examine this in some detail 
below, c) Given robust domains the population densities between each domain may well depend on the history of the 
system. The macroscopic history can be expected to condition the ensemble's set of initial microstates {1^(0)} from 
which the macroscopic material is formed. This in turn can be expected to condition the set of nonergodic domains 
{D a } that characterize the ensemble. For a macroscopic sample spanning a single ergodic domain D a , the free energy 
Gu then satisfies only a local extrema principle and thus looses much of its thermodynamic meaning. 



III. THEORY AND METHODS 



A. Equations of Motion 

We use the constant pressure Nose-Hoover equations of motion by combining the Nose-Hoover feedback mechanism 
with the so-called SLLOD or DOLLS equations of motion, which are equivalent for dilation. It is known that these 
equations of motion do not produce artifacts in the systems linear response to an external field and that to leading 
order the effect on the dynamical correlation functions is at most 0(1/N), where N is the number of particles^. The 
equations of motion are, 
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where q.; is the position, Pi is the momentum and Fj is the force on the i th particle, m is the particle mass, Ty 
is the barostat time constant, tt is the thermostat time constant, T is the input temperature, Pq is the input 
(thermodynamic) pressure and the instantaneous (mechanical) pressure is given by P(t) — Q^iLi Pi 'Pi/ m + Y^i=i F% • 
q^/ZV. Because these equations of motion have additional dynamical variables the extended phase space vector is 
r' = (r, V, ay, cut)- In order to obtain the equilibrium distribution function we first define the Hamiltonian, in 
the absence of any external fields, dilation a>v(t) = or thermostats arif) = 0, as Hq = $ + 5X^=1?^ ' Pil m i 
where $ is the total inter-particle potential energy. To proceed further we identify the extended Hamiltonian as 
He = H + ^Na T T T ksT + ^NayTyksT and then obtain the phase space compression factor 

iV 3 i~i • N 3 r\ • '"itV 

A = v • r' = V V + > > -pi + — + + — - 

= (3(H E +P V), (5) 

where the index 7 sums over the components of the Cartesian position and momentum vectors. Using the Heisenberg 
streaming representation (rather than the more usual Schrodinger representation Eq. [3} of the Liouville equation 

^ln[/(r'(t),t)] =-A(r'(t)), (6) 

we can obtain the particular time independent solution for the distribution function 

/(r') cx exp(-/3/ ) exp(~N(a T r| + a 2 v r v )), (7) 
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where Io(t) = Ho(t) + PoV(t) is the instantaneous enthalpy. The second exponential on the RHS of Eq. [7] with 
ay and qt in the argument, which has no dependence on the input temperature T or the input pressure Pq, is 
statistically independent from the rest of the distribution function, which is the standard equilibrium isothermal 
isobaric distribution. We can normalize Eq. [7|by integrating over all space to obtain the thermodynamic equilibrium 
distribution function, 

/(r') = \n T -^ exp(~N(a 2 T 4 + a 2 v 4))f Q (T, V), (8) 

Z 7T Z 

where the standard isothermal isobaric distribution function is 

f(TV)= exp(-f3(H + P V)) 

M ' ' I™dVj D drexp(-P(H +P V)y W 

It should be emphasized that the derivation of Eq. [7] says nothing about the existence or otherwise of any domains. 
These must be considered when normalizing Eq. [7] and thus Eqs. [8] & [9] are only valid in thermodynamic equilibrium. 
If we wish to use Eq. [7] outside thermodynamic equilibrium we must consider domains. 

We can also use the so called SLLOD equations of motion^ to apply strain rate controlled Couette flow (planar 
shear) to our equations of motion. The necessary modifications to the first two lines of Eq. [4] result in 

P. 

q t = 1- a v q.i + 17 lyi 

m 

pi = Fj - avPi - a T Pi - ryPyh (10) 
where 7 is the strain rate and the last three lines of Eq. [4] remain unchanged. 



B. Equilibrium Statistical Mechanics in a Single Domain 

As we have stated in the introduction, the full phase space domain includes phase points from many different 
thermodynamic phases (gases, liquids and crystals). In the thermodynamic limit this does not cause problems. To 
understand this suppose we can label microstates to be in either of two possible thermodynamic phases 1 or 2 bound by 
two phase space domains D\ and D 2 . By assumption we are not presently considering the possibility of co-existence. 
The system is assumed to be ergodic: atoms in one thermodynamic phase can in time, transform into the other phase. 
Assume the two thermodynamic phases have different free energies: G\ is the Gibbs free energy of the first phase and 
G2 is that of the second phase. For a sufficiently large N the free energy Eq. [T]is an extensive variable. We may thus 
express the partition function as the sum of contributions from the two phases 

A = + 

= e -PN 3l +e -/W» (11 ) 

where the lower case g on the second line is used to represent the intensive free energies which do not change with 
system size N in the thermodynamic limit. If g\ is less than g 2 then in the thermodynamic limit, N — > 00, the 
only significant contribution to the partition function A will be due to the "equilibrium" phase namely phase 1. Thus 
although the free energy defined in Eq. [21 is given by an integral over all of phase space D, in the thermodynamic limit 
this integral can be approximated to arbitrary precision, as an integral over the domain that includes the most stable 
phase. Suppose D\ includes only crystalline phases and D 2 includes only amorphous phases and further suppose a 
particular crystalline phase has a lower free energy that any amorphous phase. According to Eq. El we should calculate 
the free energy by integrating over all crystalline and all amorphous phases. In practice in the thermodynamic limit 
we can compute the free energy to arbitrary accuracy by integrating Eq. [21 only over that part of phase space within 
which the thermodynamically stable state resides. 

If we consider a nonergodic system that according to different preparative protocols can be formed in either phase 
1 or phase 2. After preparation, because the system is nonergodic both phases are kinetically stable indefinitely. By 
restricting the phase space integrals for the free energy to those domains that contain the kinetically stable phase 
we can compute the free energy of that phase. However, although it may be possible to formally assign free energies 
to nonergodic systems, these free energies clearly fail to satisfy any global extremum principle. As we will show 
these partition functions can be used formally to yield first and second order thermodynamic quantities by numerical 
differentiation. The metastable domain is a subset of the thermodynamic equilibrium domain which contains all 
possible atom positions including ones belonging to the metastable phase. 
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Within a single domain the system is, by construction, ergodic. Thus for almost all microstates T' a (0)£D a , 
ensemble averages, of some variable B, (B), equal time averages B, for phase space trajectories that start at time 
zero, 



Microscopic expressions for mechanical variables may be used as a test of ergodicity in nondissipative systems which 
are out of equilibrium. (Note a nondissipative system does not on average exchange heat with any thermal reservoir 
with which it has been in contact for a long time.) In the case of metastable fluids or allotropes we may introduce a 
single restricted domain and by construction the system remains ergodic within this domain. 

A gedanken experiment can be used to justify the Boltzmann weighting and the applicability of the Zeroth Law of 
Thermodynamics for such systems. Consider a double well potential with an inner and outer potential well. If the 
barrier between the inner and outer wells is much greater than fcsT, so that over the duration of observation (which is 
much greater than any relaxation time in the ergodically restricted subsystem) no particles cross the barrier, then the 
system considered as a double well system, will be, by construction, nonergodic. For systems composed of particles 
that are solely found in the inner potential well, our hypotheses are that the distribution of states in the inner well 
will be given by a Boltzmann distribution taken over the inner domain only and that if such a system is in thermal 
contact with a body in true thermodynamic equilibrium, then the temperature of the ergodically restricted system 
must equal that of the system in true thermodynamic equilibrium. We can justify these hypotheses by considering a 
fictitious system that only has the inner potential well and in which the potential function is positive infinity for all 
separations that are greater than the inner well (this includes the position of the outer well) . In accord with Gibbsian 
statistical mechanics the distribution of states is canonical over this (single well) potential. Furthermore the Zeroth 
Law of thermodynamics will apply to this single well system. Now if we dynamically generate the outer well, all the 
particles locked inside the inner well cannot "know" that the outer well has been formed so their dynamics will be 
completely unchanged by the time dependent generation of the new outer well. The generation of an inaccessible 
outer well will not alter the distribution of states in the inner well nor will it cause any flow of heat to the equilibrium 
heat bath surrounding the system. This provides a compelling physical justification for our domain hypotheses over 
a single ergodic sub domain of phase space. 

In order to recover many of the basic relationships of Gibbsian statistical mechanics it is also necessary that the 
system appears to be in dynamical equilibrium, i.e. | f(F a ,t) — f(T a ,t + r ) |< e,V T a . G D a , for some small e, over 
the longest observation time r c . We use the definition of the partition function Eq. [2] as before but now the domain 
D a in the integral is over a single contiguous hypervolume in the configuration space of the generalized position 
coordinate q and volume V. The domain over the generalized momentum p and multipliers ay and ar remains 
unchanged. We then obtain the Gibbs free energy by use of Eqs. [U&El Thus far all we have altered is our definition 
of the domain. In changing the definition of the domain we have opened a potential problem for Gibbsian statistical 
mechanics. If we change the temperature or the pressure of the system the domain may also change. If the domain 
changes this may make a contribution to the derivatives of the partition function, Eq. [21 and the direct connection 
with the standard outcomes of macroscopic thermodynamics will be lost. Thus the domains need to be robust with 
respect to changes in thermodynamic state variables. 

There are three means by which a system could have robust domains. The first and most obvious is that the domain 
does not change when the pressure or the temperature changes, 8D a (X,Y)/dX = 0, where X is a thermodynamic 
state variable and Y is the other thermodynamic state variables. When we lower the temperature only the inverse 
temperature j3 in Eq. [7] changes and when we change the pressure only the parameter Pq changes. If the domain's 
boundary is determined by a surface on which Io(T') always has a very high value it will remain unchanged under 
infinitesimal changes in Pq or (3. We will refer to a surface domain that doesn't change with the state variables as 
completely robust. The second way is that the distribution function is always identically zero on the domain boundary, 



where S a is the surface of the domain D a . Because the domain is contiguous (required for it to be ergodic) it must 
have a single connected surface. Such a domain will be robust. The third way the domain can be robust is less 
restrictive and allows for the possibility that the domain does change when the thermodynamic variables are changed 
substantially. If SX is an infinitesimal change in a thermodynamic state variable then, 




(12) 



fo(T,V) =0 VfeS, 



(13) 



6D a (X + SX, Y) = 6D a (X, Y) + 0(SX) 



n 



(14) 
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where Y denotes the other thermodynamic state variables, we require that n > 2 for first order thermodynamic 
property formulae to be correct n > 3 for second order property formulae to be correct etc. Obviously this third way 
will be satisfied in the first two cases as well. 

Later in the paper we will introduce an independent test of domain robustness. However, if a system was not 
robust then we would expect that small changes in the state variables would change the macroscopic properties of 
the sample permanently - it would be as though the preparation history of the sample was continuing even for small 
changes in the state variable. Quite obviously if we produce huge changes in the state variables we will of course 
permanently change the properties of the system because we permanently deform the ergodic domain. Experience 
shows however that very many nondissipative nonequilibrium systems are quite robust with respect to small changes in 
state variables. All that is required for fluctuation formulae for first second and third order thermodynamic quantities 
to be valid, is that the domains be unchanged, to first second or third order, by infinitesimal changes in the state 
variables. Obviously a robust domain is an ideal construct. However on the typical time scale of interest, which is 
usually orders of magnitude less than the time scale on which the system will change to a new phase of lower free 
energy, this can be a very good approximation. 

We are now able to recover most of the standard results of Gibbsian equilibrium statistical mechanics. For example 
we may calculate the enthalpy (I) a from the partition function A, Eq. as 

_ J J Da dV dTI (T 1 V)exp(-f3I (T,V)) 

jj Da dVdTexp(-0I o (T,V)) ~ a - (15) 

Here we are considering an ensemble of systems which occupy a single ergodic domain D a . Since this domain is self 
ergodic the ensemble average is equal to the corresponding time average. 

The term on the RHS of the second line is obviously the average value of the instantaneous enthalpy, io(r, V) = 
Hq(T,V) + PqV, obtained by using the equilibrium distribution function, Eq. [9] or equivalently Eq. [8] with the 
integration limits restricted to the domain D a , where Po is the externally set thermodynamic pressure. We can also 
obtain expressions for the average volume (V) and the constant pressure specific heat cp by taking the appropriate 
derivatives of the partition function Eq. [2 In other ensembles we can use the same procedure to find other variables 
e.g. the internal energy, the average pressure and the constant volume specific heat in the canonical (N, V, T) ensemble. 

An important outcome is that this description remains compatible with macroscopic thermodynamics. Here the 
Gibbs free energy is defined as, 

G = U - TS + P (V) , (16) 

where U = (Hq) is the internal energy and S is the entropy. If we take the derivative of Eq. [16] with respect to one 
of the isobaric isothermal ensembles conjugate variables (N, Po, T) while keeping the others fixed we obtain 

\dPo) Nt 

We now write down the microscopic equilibrium equation for the Gibbs entropy 

S = ~ kB J L dVdT M T ^ ln M^,V). (18) 

It is an easy matter to show that Eqs. [J, [2] & [18] are consistent with the two derivatives given in Eq. [171 Given our 
condition of a robust boundary we thus have a form of Gibbsian statistical mechanics for metastable states which 
remains in agreement with macroscopic thermodynamics. 

C. Multiple Domains and Nonergodicity 

We now wish to consider an ensemble of systems which is prepared from an initial ergodic (usually high temperature) 
equilibrium ensemble. There is some protocol V, which involves a temperature quench or a sharp pressure increase 
etc, which breaks the ensemble into a set of sub-ensembles characterized by different macroscopic properties. After 
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(V) 



(17) 
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the protocol V, has been executed we allow all the ensemble members to relax to states which are macroscopically 
time independent - to within experimental tolerances. We assume that the ensemble can be classified into a set 
of sub-ensembles {a, a = l,Nu} whose macroscopic properties take on Njj distinct sets of values. For the longest 
observation times available a macroscopic system classified as an a system is not observed to transform into a 
system, and vice versa. The full ensemble of systems is thus non-ergodic. However, in each individual sub-ensemble, 
say sub-ensemble a, the constituent members are ergodic (by construction). Thus we can partition the full phase 
space into a set of domains, {D a }. 

From the arguments given above (in section B), after the relaxation of initial transients, we expect to observe a 
Boltzmann distribution of states within an individual domain which is therefore independent of the quench protocol. 
However the distribution between domains cannot be expected to be Boltzmann distributed and will instead be 
dependent on the quench protocol. Within a given ensemble the proportion of ensemble members ultimately found 
in domain D a is given by a weight w a (V) which is subject to the constraint 



Nd 



a=l 



We can calculate the full ensemble average of some macroscopic property B as, 

/m ^ I Id dVdTB(T,V)exp(-/3I ) 

a=l 



J J D dVdrexp{-pI ) 



(19) 



(20) 



Since the full ensemble of states is nonergodic the phase space breaks up into disjoint domains which in themselves 
are ergodic. Thus each domain may be identified by any point in phase space, (r, V), that is a member of it so the 
subscript a is a function of the phase vector a(T,V) allowing the following expression for the distribution function 



Nr. 



f(T,V) = J2^s(T,D a )f a (T,V) 



(21) 



where s(T, D a ) = 1 if TeD a and s(T, D a ) = otherwise, and 

exp(-/?/ (r,V)) 



fa(T,V) 



J™dVf D drexp(-/3/ (r,V))' 



(22) 



The entropy is given by S = —k B dV J D dTfln(f) and using Eq. [21] we have the following expressions for the 
multidomain entropy, 



Nd r „oo r ~i N d N D 

S = -k B ^ w a / dV dr/c, ln(/ Q ) +ln(w Q ) = ^ w a S a - k B V] w a \n{w a ) 
a =i Lio Jd c , J a=1 a=1 



(23) 



The term — 5^a=i ^BW a ln^Q.) = Sd is the inter-domain entropy, which is maximized by an even distribution of 
ensemble members over all domains, while S a is the intra-entropy of domain a considered as a single TV-particle 
system. 

If we substitute Eq. [22] into Eqs. [23] we find that, 



S - S D +T- 1 Y J w a {I ) a + k B Y J wJn f dV f drexp(-/3J (r, V)) 

a =i a=i J o JD a 

Nd poo /> 

= Sd+T- 1 (I) +k B Y]wJn I dV drex P (-/?/ o (r,V0), 

a=l JO JD a 



(24) 



where (B) a = J J D dV dT B(T)f a (T,V). Combining Eq. [16] with Eq. [24]we obtain the following expression for the 
Gibbs free energy 



Nr 



G = -k B Tj2 



lu/ / dVdTexp(-(3I ) -\n(w a ) 



N D 

w a G a - S D T. 

a = l 



(25) 
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It is easy to verify that if we hold the local domain weights fixed and then vary the temperature or the pressure 
that Eqs. [23] & [25] are compatible with Eqs. [TT] This means that if we have a fixed number of robust domains, 
whose population levels or weights are non-Boltzmann distributed, Eqs. [25] & [23] provide a direct microscopic link to 
standard macroscopic thermodynamics. On the extremely long time scale the weighting functions w a may vary and 
the system will tend towards the direction where the free energy Eq. [25] is reduced. Without the inter-domain entropy 
term, So, Eq. [251 would be minimized when the domain with the lowest free energy has all the ensemble members in 
it. It turns out that Eq. [25] is minimized when all the domain weights are Boltzmann distributed, ie when 

JJ Da dVdrexpj-pio) 

W a = K? ■ (26) 

E^Jf Dff dVdTexp(-pi ) 

Here (i.e. upon obeying Eq. [26]) the entropy and free energy given by Eqs. [23] & [25] coincide with the standard 
equilibrium expressions so the free energy must be a minimum. To prove this we use Eq. [25] and we remove the first 
weight wi = 1 — J2a=2 Wa ' so that the constraint Eq. [19] is respected while the remaining weights are independent. 

This means that the free energy can be written as G(l — w a , W2, W3, ...wn d )- The constrained partial derivatives 



are then 

dG 



dw D 



a=2 



dG dwi dG dG dG 

a > 2 



dw\ dw a dw a dw\ dw a 
dG dG 

a>2, 



dwi dw a 

= -Gi - kTln( Wl ) - k B T + G a + k B Tln(w a ) + k B T, a>2. (27) 
Using the fact that for Boltzmann weights, Eq. (26] 

w a = exp[/3(G eq - G a )\ (28) 
where G eq = — k B Tln X) a =o / Id dV dT exp(—f3Io) is the equilibrium free energy, we find that at equilibrium, 

dC 

' ' 1 = 0, a > 2. (29) 



dw, 

It remains to prove that this is indeed a minimum. Using the same approach for treating the constraint, we continue 
making the first weight a function of all the others, and obtain 



d 2 G 



dw a dw 7 



d 2 G d 2 G d 2 G d 2 G , „ 

(30) 



dw\ dwjdwa dw 1 dw\ dw\dw a 
Using the Boltzmann weights it is easy to show that 

d 2 G kT 

-—2 = — = k B Texp[(3{Gi - G eq )] 
dwf w\ 

d 2 G d 2 G 



duo \d w a dw-fdwi 
d 2 G „ k B T 







dwjdw a w, 



Sja = 6 ia k B Texp[(3(G a - G eq )]. (31) 



G 



is positive definite, thus we have proved the free 



From these results it is easy to see that the Hessian matrix 75 

energy to be a minimum for the case of equilibrium. 

We can use a knowledge of the multiple domain thermodynamic potential, Eq. [25] to compute averages. As an 
example we consider the average enthalpy again, 



, M 2 d(3G 

(J ) = -k B T 



dT 



(32) 



,Pn,N 



Eq. [32] can easily be derived from Eq. [16] and is in essence the same as the first line of Eq. [HI It is straightforward 
to see that upon using Eq. [25] to calculate the derivative in Eq. [32] one obtains the average enthalpy as given by Eq. 
[20l One can do the same for other quantities such as the specific heat etc. 
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D. Application to Molecular Dynamics Simulation 



We are now in a position to test various outcomes, using computer simulation, which may be derived by drawing 
on the previous material. If we start with an ensemble of systems which are initially in equilibrium at temperature 
T = Xb and then at time t = we quench them by setting T = T\ we can solve the Liouville equation Eq. [6] to obtain 

f(T(t),a v (t),a T (t),t) = /(r(0W(0),« T (0),0) 

x exp^/oW-ioW)] (33) 

where /3i = X/ksTi. This nonequilibrium distribution function, valid for t > 0, explicitly requires the solution of the 
equations of motion and is of limited utility. However it allows the identification of a formal condition to identify 
the amount of time, which must elapse after the quench, before Eqs. [TJ [H [20] or [21] can be applied to the ensemble. 
That is the quantity Io(t) — H (t) + PoV(t) must be statistically independent of Iq(0). Thus we are interested in the 
correlation function 

Ci m)MQ)) = (/oWJo(0))-(/oW)(/o(0)) (34) 

wo 



where 



Ci,o = ^((/oW 2 ) - {Io(t)f) ((/o(0) 2 ) - </o(0)) 2 ). 



(35) 



This function will equal 1 for a perfectly correlated system, -1 for a perfectly anticorrelated system and for an 
uncorrelated system. When we consider an ensemble of systems, occupying the various domains to different degrees, 
we see that Eq. [34] may not decay to zero given the trajectories are unable to leave their domains. If the transients, 
due to the quench, fully decay Eq. [34] will decay to zero and Eqs. [7] & [9] will become valid for the ensemble. 

If the correlation function, Eq. [34] decays to a plateau then we may have a situation where Eqs. l20l and I2T1 are 
valid. It may be that the material can still slowly age, due to processes, that occur on a time scale which is longer 
than the one we are monitoring. For a glass we expect that correlation function Eq. [34] will not fully decay on a 
reasonable time scale. If we give the system time to age, such that it appears to be time translationaly invariant, and 
then compute the following correlation function, 

C 2 (/ (r),/o(0)) = Ztoilomm -(Io(t)) a (Io(0)) a (M) 



C 2 



where 



C2,o — 



E 

a=0 



\J ((Io(t) 2 ) a (Io(t))l) (Vo(0) 2 ) Q - (l (0))l) , (37) 



we may observe a full decay. If this occurs Eqs. [20] and [21] will be valid. To compute this correlation function 
N t trajectories are produced and for each of these the averages, (B) a , where B is an arbitrary dynamical variable, 
appearing in Eq. [36]are approximately obtained by time averaging. When the system is ergodic and time translationaly 
invariant these two correlation functions, Eqs. [34]&[36l will give the same result. However for a nonergodic system C2 
will decay to zero on a reasonable time scale while C\ will not. Rather C\ may decay to some plateau on a reasonable 
time scale and then decay on a much slower time scale. The preceding sections then rest on this clear separation of 
time scales in the correlation function Eq. [34j For metastable fluids and allotropes this separation is so extreme that 
we probably cannot observe, even the early stages of, the later slow decay on any reasonable experimental time scale. 
Further for these systems there will be only a single domain and thus they appear ergodic. For glasses some signs of 
the later decay can often be observed, however it is still very much slower than the initial decay. In the field of glassy 
dynamics the initial decay is often called the (3 decay and the slower long time decay is often called the a decay. As 
the glass is further aged this second stage decay is observed to slow down dramatically while the early decay does not 
change very mucbM. 

To allow the hypothesis of local phase space domains to be tested we will now introduce several relations whose 
derivation draws upon the equilibrium distribution function Eq. [9] We will also discuss the effect of the phase space 
domains on these relations. 

First we introduce the configurational temperature 



(38) 
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where F is a 3iV dimensional vector representing the inter-particle forces on each atom F = — V$. This relation is 
easily derived from Eq. [9] (set B to V • F, integrate by parts and drop the boundary terms) and will remain valid for 
nonergodic systems where Eq. [26] is not obeyed. The relation involves the spatial derivative of the force, which is not 
zero at the cutoff radius for the potential we use in our simulations (see below) . This along with finite size effects can 
result in a small disagreement between Eq. [38] and the kinetic temperature for our system when in equilibrium. 

Equilibrium fluctuation formulae may be easily derived from Eq. [2Ql see ref.— £ for some examples of how this is 
done. We consider how the average enthalpy changes with the temperature at constant pressure, 



Cp = (V) c Po 



d_ 

dT 



N,P 



(lo) = | {(l 2 o) (hY 



(39) 



where cp is the constant pressure specific heat. The calculation of the RHS of Eq. [39] using the ensemble average 
given by Eq. [20] is subtle. If we simply calculate (Iq) and (Iq) with the use of Eq. [20] and then plug the results into 
Eq. [39] we obtain, what we will refer to as, a single domain average which does not give us the correct change in 
the average enthalpy for a history dependent equilibrium ensemble. This is because the average (7q) superimposes 
across the different domains while the quantity (Iq) 2 contains spurious cross terms. If we derive the heat capacity by 



taking the second derivative of the Gibbs free energy for the multiple domain distribution function, 
with respect to the temperature, 



Eqs. 



[&l 



r 9 



k B T< 



d(5G 



N,P ,w a 



dT 



(40) 



N,Pn,W a 



we obtain the following 



(3 Nd 

Cp = (V) c Po = — w 

Q = l 



(41) 



where the quantity (Jq) — (Io) a ls obtained for each domain separately (here (. . .) a represents an average taken where 
all ensemble members are in the a th domain). We will refer to this as a multidomain average which is consistent 
with the nonergodic statistical mechanics and thermodynamics that we have introduced here. It is obvious that 
both a single and multiple domain average will give the same result in the case of thermodynamic equilibrium and 
metastable equilibrium (single domain). The transition from the single domain average producing the correct result 
to an anomalous result is symptomatic of an ergodic to a history dependent nonergodic transition. If we consider 
a large macroscopic system (the super system) to be made of N s independent subsystems the multidomain average 
remains self-consistent. To see this we apply Eq. [41] to fluctuations in the super system and then we inquire how 
this relates to fluctuations in the subsystem. The enthalpy of one instance of the super system will be given by 
I s = J2a=i I®- I n principle an ensemble of super systems can be prepared by applying the same history dependent 
macroscopic protocol to all members of this ensemble. Due to the statistical independence of the subsystems, upon 



taking an ensemble average of super systems, we have (I a Ii 



Pis 



(I a )s (Ip)s f° r all a 7^ /3. Here the average 



is 



taken over the ensemble of super systems and the a th subsystem in each super system is identified by its location. It is 
then easy to show that the specific heat obtained from the ensemble average Eq. [41] of the super system is equivalent 



to that obtained from the subsystem due to the two quantities (lg) s = \ (Yl a =i ^ a J ) anc ^ (I s ) s = \l2 a =i ^ a 

possessing identical cross terms which cancel each other out (as a result of the independence of the subsystems) upon 
applying Eql4Tl 

If we ignore finite size effects, due to assuming the equivalence of ensembles, the constant volume specific heat is 
related to the constant pressure specific heat by the equation 



Cv = (V) c v = Cp 



dV 




(OH 




1 dV 




dV 


dT 


T 


V ap 


t/ 




r) 


dT 



(42) 



We may also obtain an expression for the constant volume specific heat cv by deriving equilibrium fluctuation formula 
for each of the derivatives appearing in Eq. [42]in lieu of directly measuring them. We may then obtain a single domain 
expression for cy which does not work for the history dependent glass and also a correctly weighted ensemble average 
(multidomain average) which does. This is completely analogous to what has been shown in detail for cp. As the 
equations are unwieldy, and their derivation (given an understanding of the cp case) is straightforward, we will not 
reproduce them here. 



11 



E. Test of Domain Robustness: Transient Fluctuation Theorem 



The application of the Evans-Searles Transient Fluctuation Theorem ^ 15 ' 16 to the systems treated in this paper 
provides a sharp test of the assumptions used to develop the theory given in this paper. The Theorem describes 
a time reversal symmetry satisfied by a generalized entropy production, namely the so-called dissipation function. 
The precise mathematical definition of this function requires a knowledge of the dynamics and also of the initial 
distribution function. The three necessary and sufficient conditions for the Fluctuation theorem to be valid are that 
the initial distribution is known (here we assume the distribution is Boltzmann weighted over some initial domain of 
phase space) , that the dynamics is time reversible (all the equations of motion used here are time reversible) and lastly 
that the system satisfies the condition known as ergodic consistency. When applied to the systems studied here this 
requires that the phase space domains should be robust with respect to the sudden changes imposed on the system 
and that the number of inter-domain transitions remain negligible on the time over which the theory is applied. If 
any one of these three conditions fails then the Theorem cannot be applied to the system and the corresponding 
fluctuation relation will not be satisfied^. 

We can use the fluctuation theorem to obtain relations for how the system responds upon suddenly changing the 
input temperature or pressure for a system, which is initially in equilibrium as specified by Eq. [7] or [8] Firstly we 
consider a change in the pressure, while holding the temperature fixed, by changing the input variable Pq in Eq. [4] 
(thermodynamic pressure) to Pq — P2 at time t — for a system initially in equilibrium with Pq — Pi- The probability 
density p(AV = A) of observing a change in volume of AV(t) = V(t) — V(0) relative to a change of equal magnitude 
but opposite sign is then given by 

p^nt)=A) _ exp{p{P2 _ Pi)Ay m 



p(AV(t) = -A) 

To derive this expression we have had to assume that the intra-domain populations are Boltzmann distributed accord- 
ing to Eq. Ergodic consistency requires that for any initial phase space point T(0) that can be initially observed 
with nonzero probability, there is a nonzero probability of initially observing the time reversal map M T of the end 
point r(i), (ie VT(0) such that /(r(0),0) ^ 0, f(M T (T(t)), 0) ^ 0). This condition obviously requires that the phase 
space domains remain robust and the number of inter-domain transitions remain negligible for at least a time t, after 
the pressure (or temperature) quench. 

If we sample all or our initial t = 0, Po = Pi configurations from the one trajectory which remains locked in a 
single domain even after the quench, we expect Eq. [43] to be valid. If we prepare an ensemble of initial configurations 
using the same protocol we still expect Eq. [43] to remain valid even with different domain weightings Wi, as defined 
in Eqs. [2Q]&[2lJ provided the domains are robust over the time t appearing in Eq. [43] Note that suddenly reducing 
the pressure by a very large amount could result in a breakdown of the robustness condition. Eq. [43] may be partially 
summed to obtain what is referred to as the integrated fluctuation theorem 

^|^ = (exp(/3(P 2 -P 1 )A^)) AV<0 . (44) 

For the case where we change the input temperature T in Eq. [4] while holding the input pressure Po fixed we obtain 
a relation for fluctuations in the extended instantaneous enthalpy = Hsit) + PoV(t). We start with a system 

initially in equilibrium at temperature T — l/(fcs/3i) and we then subject it to a temperature quench by changing the 
input temperature in Eq. [4] to T = l/ftsfa), at time t — 0, while holding the input pressure fixed. The probability 
density p(Ais(t) — A) of observing a change in instantaneous enthalpy AJe(*) = /#(£) — Ie(0) relative to a change 
of equal magnitude but opposite sign is then given by 

P (AI E (t)=A) _ expm _ my (45) 



p{AI E {t) = -A) 



Note that if we suddenly increase the temperature by a very large amount we could expect to violate the robustness 
or the negligible inter-domain transition condition. In common with Eq. [43] we expect that this expression will be 
valid when all initial configurations are sampled from a single common domain and also when sampled from multiple 
arbitrarily populated domains under the assumption that the domains are robust and the number of transitions are 
negligible over time t. This equation may also be partially summed to obtain 



(46) 
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Figure 1: A logarithmic plot of the self diffusion coefficient for both species A and species B particles as a function of the 
separation parameter T — T g with T g = 0.435. 



IV. SIMULATION DETAILS 



For our simulations we use a variation on the Kob Andersen glass former™ featuring a purely repulsive potential^. 
The pairwise additive potential is 



Uij(nj) = 4e Q/3 



fe) 12 -fe) 6 + i 



v m < ^2 



iHj(rij)=0 V rij > \/2a a /3, (47) 

where the species identities of particles i and j, either A or B, are denoted by the subscripts a and (3. The energy 
parameters are set ebb — 0.5eAA, £ab = 1.5 and the particle interactio n distances <jbb = 0.88ctaa, oab = 
0.8 a aa- The energy unit is €aa, the length unit is ctaa and the time unit is y/m&AA/ £ aa with both species having 
the same mass m. The composition is set at X = Nb/Na = 0.2, the number of particles are N = Na + Nb — 108, 
the pressure is set to Pq — 14 eaa/c 3 and the temperature unit is €aa/^b- The time constants are set at ry = 
and tt = VN, Note that the energy parameters are slightly different to the potential we used in refpi^. The equations 
of motion were integrated using a fourth order Runge-Kutta method^. The time step used was dt = 0.002 and 
sometimes dt = 0.004 for very low temperatures. 

From previous work on binary mixtures we know the basic reason why this system is vary reluctant to 
crystalliz o 20 i 21 i 22 i 23 . The chosen nonadditivity of the species A-B interaction makes the mixture extremely misci- 
ble; consider the present value of gab = 0.8<taa relative to the additive value of oab = 0.94 ctaa- This effect 
dominates over the choice of the energy parameters. Due to this extreme miscibility the relatively large composition 
fluctuations necessary, about the average composition of X = 0.2, to form the crystal phases (either the pure species 
A, X — 0, FCC crystal or the binary, X — 0.5, CsCl crystal) are strongly suppressed and crystallization is strongly 
frustrated. 

We identify the nominal glass transition by calculating the diffusion coefficient as a function of temperature. This 
is shown for both species in Fig. Q] on a logarithmic plot demonstrating how the diffusion coefficient approaches 
zero critically, D a (T) oc (T — T g ) b where b is the critical exponent, with a nominal glass transition temperature of 
T g = 0.435. It would be, perhaps, more customary to obtain a nominal glass transition temperature by analyzing 
the critical divergence of the viscosity. Given that the Stokes Einstein relation is strongly violated upon approaching 
the glass transition one might be concerned that this would give a very different result. However the violation of 
the Stokes Einstein relation can largely be attributed to the exponent b being different between the viscosity and the 
diffusion coefficient rather than the nominal glass transition temperature T^2£. 
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Figure 2: a) The instantaneous enthalpy correlation functions as defined by Eqs. [34] & [36] as a function of logarithmic time for 
the temperature T = 1. The calculation of the function was started after the fluid had been given time to equilibrate. The 
strong agreement between the two correlation functions is indicative of ergodicity. 

b) The instantaneous enthalpy correlation functions as defined by Eqs. [34] & [36] as a function of logarithmic time for the 
temperature T = 0.4. The calculation of the correlation function was started at various times after the quench, all approximately 
at t = 8 x 10 5 , in an attempt to age the system. The difference between the correlation functions is indicative of nonergodicity. 
Notice that C2 reaches a near full decay between t = 1 & 10, while Ci reaches a non-decaying plateau. 



V. RESULTS AND DISCUSSION 



The correlation functions given in Eqs. [34] & [36] were calculated from ensembles of 100 independent simulations at 
the two temperatures given in Fig. [2](T = 1 and T = 0.4). In all cases the systems were subject to an instantaneous 
quench, from an initial equilibrium at T = 5, by changing the value of the input temperature in Eq. [4] The system 
was then run for a significant time, in the case of the glass ensemble T age > 8 x 10 5 , in an attempt to age it. Of course 
the longest time that can be accessed in a molecular dynamics simulation is rather short, and so the system is not 
very well aged, but we are still able to meaningfully treat it as a time invariant state. Each of the 100 independent 
simulations was interpreted as being stuck in its own domain D a and the correlation functions were calculated for 
each of these domains using time averaging. The time averaging was approximately 100 times longer than the longest 
time t = 800 that the correlation functions were calculated out to. Obviously in the limit of an infinite number of 
independent simulations and the case where the domains are robust we will obtain the exact multidomain average 
given by Eq. [2Ql We assume our limited ensemble of simulations is representative of this. The data from each domain 
(independent simulation) was then used to obtain the correlation functions Eqs[34] & [36] as seen in Fig. EJ At the 
higher temperature T = 1 it can be seen that the two functions are equivalent demonstrating how the system is 
ergodic. It can also be seen that the correlation function has decayed on a time scale of t ~ 10 which is therefore (by 
Eq. [33]) the time scale on which the ensemble becomes accurately represented by Eq. [9] with only one domain D which 
does not necessarily extend over all phase space. The oscillations, which can be seen in the correlation function, are 
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due to both ringing in the feedback mechanisms of Eq. [4] and the frequency dependent storage component of the bulk 
viscosity. The statistical uncertainty in the correlation function becomes larger than these oscillations somewhere 
between a time of t — 1 & 10. If we constructed an experiment where the pressure was regulated by a piston and a 
spring, the correlation functions, Eqs[34]&[36] would depend on the details of the piston and spring parameters in a 
similar way to the simulations dependence on the details of the feed back mechanism. 

When the system is quenched to the lower temperatures (T = 0.4) ergodicity is lost and we obtain a glass. The 
complete decay of the first correlation function, C\, Eq. [34] may no longer occur because the individual trajectories 
remain stuck in local domains which have different average values for (Io) a - This is similar to the much studied 
density correlation function^ (the intermediate scattering function) which decays to a finite plateau for a glass or 
more generally a solid material. On the other hand there is nothing to stop the second correlation function, C 2 , Eq. 
l36l from fully decaying when the system is nonergodic. If the second correlation function fully decays whilst the first 
is only able to decay to a plateau then we have a situation where Eqs. [21] & [20] are valid as can be seen from Eq. [33] In 
Fig. [2] it can be seen that C\ does indeed fail to decay while Ci comes very close to fully decaying at the time where C\ 
reaches the plateau. The reason C2 doesn't fully decay here is due to the fact that the inter-domain transition rates, 
while small, are not exactly zero. If the system had been aged more extensively this problem would be significantly 
reduced. This effect is exacerbated by the time averaging, used to form the averages for each trajectory, being two 
orders of magnitude longer than the longest time the correlation function was calculated to. The effect of the state 
slowly evolving due to finite inter-domain transition rates is too small to seriously compromise the modeling of the 
system as obeying Eqs. [20] & EH and thus we have obtained direct evidence for the validity of these equations. The 
height of the plateau for C\ will depend on the history of the system, i.e. the protocol used to prepare the ensemble. 

We move on to a comparison between the kinetic temperature and the configurational temperature, the results of 
which may be seen in Fig. [3^)- The input temperature ranges from T = 3 to T = 0.3. Also shown are results for the 
system, undergoing constant planar shear, Eq. [TO] with a strain rate of 7 = 0.5. At the higher temperatures we see 
a very small relative discrepancy between the two types of temperatures, which we attribute to the discontinuity in 
the first spatial derivative of the inter-particle force at the cutoff radius and to finite size effects. These effects appear 
to diminish a little at lower temperatures. For temperatures above T = 1.5 the chosen strain rate has no significant 
effect on the configurational temperature indicating that our system is in the linear response domain^. At the lowest 
temperatures, well below the glass transition temperature, we observe good agreement between the configurational 
and input temperatures for the system without shear. This provides further evidence of our assertion that the system 
obeys Boltzmann statistics in the glass Eqs. [2Q]&[2T] However, at low temperatures, the system that is undergoing 
shear shows an increasing relative discrepancy between the two temperatures. At low temperatures the system leaves 
the linear response domain^ demonstrating the fundamental difference between the nonequilibrium distribution of 
the history dependent glassy state and that of a strongly driven steady state. If we drive the system hard enough, 
at any given temperature, we can always make a disagreement between the two types of temperature due to the 
steady state no longer being accurately represented by a Boltzmann distribution i.e. due to a break down in local 
thermodynamic equilibrium. When the system is not driven by an external field we have been unable to observe any 
difference in the two temperatures by deeply supercooling a glass forming mixture apart from the initial transient 
decay immediately following the quench, which falls off surprisingly rapidly. 

In Fig [3b) results are presented for the heat capacities (the specific heat multiplied by the volume) at both constant 
pressure Cp and constant temperature Cy . Details of the protocol used to obtain this data is given in the end 
noteSi. The estimates from the multidomain average are compared with those from the single domain average. The 
results from the multidomain averages exhibit the well-known peak, which is a signature of the onset of the glass 
transition, and has been observed directly by calorimetry in many experiments on real glass forming materials^. 
The temperature, where the peak is observed, depends on the history of the system. No peak is observed for the 
single domain averages which continue to increase as the temperature is lowered. While the two methods for forming 
averages give the same results at temperatures above the peak, they diverge at temperatures below the peak. It is 
the multidomain average that gives results consistent with the actual calorimetric behavior of the system. This may 
be seen in the figure by comparing the data which has been computed by numerically differentiating the enthalpy 
using central difference. At the peak neither the central difference (due to rapid rate of change) or the multidomain 
average data (due to a lack of domain robustness) are reliable and they show significant differences. However below 
the peak they once again show quantitative agreement providing strong evidence that the domains are robust in this 
region. If we substantially increase the duration the time averages (for each domain) are constructed on, the peak will 
be shifted to lower temperatures as previously shown2I. This requires the time average to be constructed over some 
two orders of magnitude more time than the decay time for the correlation function in Eq. [34] This is necessary in 
order to obtain enough independent samples for a meaningful estimate of the variance, of the instantaneous enthalpy 
appearing in Eq. [39] For a large macroscopic system we would expect that the specific heat measured over the entire 
ensemble would differ very little to that measured from any one of its members. We are now in a position to make 
an unambiguous interpretation of the peak in the specific heat. The peak is observed at the temperature where the 
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Figure 3: a) The relative difference between the kinetic temperature T (controlled directly by the Nose Hoover thermostat) 
and the configurational temperature given by Eq. [38] as a function of the kinetic temperature. The results for the system 
undergoing Couette flow (Shear) are for a constant strain rate of 7 = 0.5. 

b) The heat capacity calculated using the equilibrium fluctuation formula by the single domain averaging method Eq. I39I and 
the multidomain method Eq. [41] Also shown is data obtained by numerically differentiating the enthalpy by central difference. 
At the temperatures above the peak the three types of averages give very similar results. Also shown is equivalent equilibrium 
fluctuation formula data for the constant volume specific heat Eq. [42] 



system leaves metastable equilibrium and enters a history dependent state that requires averages to be computed by 
Eq. [2Q]rather than by direct use of Eq. [9] The calculation of both (Io) a and (7q) wm be different for each domain. 
If we use time averaging to calculate these variables on a time scale that falls within the plateau region for Eq. 03 
see Fig. [2b ; then the amount of time chosen to form the average is not critical. The peak occurs because the various 
ensemble members have become locked in local domains on the time scale that we are able to access. Near the peak 
itself these domains are not expected to be robust. 

The multidomain average, Eq. HTJ gives the heat capacity for a glass with robust domains. At the lowest tem- 
peratures the heat capacity reaches the beginning of a plateau, Fig. [3b. For the constant volume heat capacity Cy 
this plateau (within uncertainties due to finite size effects) has a value that is consistent with the Dulong-Petit law«2&, 
as would be expected for an amorphous solid where the potential energy surface can be modeled as harmonic upon 
transformation to the orthogonal independent basis set. This is exactly what we would expect from our local domain 
model at low temperatures. 

Testing the integrated transient fluctuation theorem (ITFT) for a sudden pressure change Eq. [44] and a sudden 
temperature change Eq. [46] provides further evidence that the Boltzmann distribution may be used to accurately 
describe intra-domain statistics in the glassy state and also that in a properly aged glass, the domains are robust 
with respect to the pressure and temperature changes studied here, Fig. [4] These equations remain valid whether we 
subject an ensemble of simulations (multidomain) to a quench or we sample from a single trajectory (single domain), 
which remains stuck in a single domain. The accuracy with which these relations are satisfied is powerful independent 
evidence for the applicability of our assumptions to the systems studied here. The fact that over the times shown in 
Fig. [H the ITFT does indeed yield correct results directly implies that, within experimental tolerance of the data, 
the phase space domains must be robust and the number of inter-domain transitions must be negligible. Unlike the 
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Figure 4: Results from applying the fluctuation theorem to the glass for a) a sudden pressure change, where the symbols are 
p(AV > 0)/p(AV < 0) and the solid line is(exp(/3(P2 — Pi)AV)) AV<0 and b) a sudden temperature change, where the symbols 
are p(AIe > 0)/p(AIe < 0) and the solid line is(exp((/32 — /3i)AIe)) ai <0 . Results from an ensemble of independent initial 
systems and in addition from a single initial trajectory (with a time of t = 5 being computed between each transient trajectory) 
are shown for a total of 10 5 pressure or temperature changes. The serial results have been shifted up, for clarity, by adding 0.2 
to the data 

specific heat fluctuation formula this requires that the domains are robust to finite changes of the state rather than 
infinitesimal changes. Thus, given that we have aged the glass sufficiently that domains are robust the number of 
transitions are negligible over the longest time the fluctuation formulae are computed, we expect Eqs. |44]&[46]to be 
correct. If we wished to apply the steady state fluctuation theorem matters become more difficult. 

VI. CONCLUSIONS 

We have presented a rigorous development of statistical mechanics and thermodynamics for nonergodic systems 
where the macroscopic properties are sensibly time independent and the phase space for the ensemble is partitioned 
into robust domains. Using computer simulation we have carried out various tests on a glassy system and shown that 
apart from the immediate vicinity of the glass transition, the computed results are consistent with our theory. While 
the intra-domain populations are individually Boltzmann distributed, the inter-domain populations are not. 

A correlation function whose decay to zero requires global Boltzmann weighting, has been derived and it has been 
shown that it decays on a reasonable time scale for ergodic systems but not for nonergodic systems. A second 
correlation function which decays to zero if the intra-domain populations are Boltzmann distributed but globally the 
inter-domain populations are not, has also been defined. We have developed expressions for obtaining averages in a 
multiple domain ensemble and shown how single domain averages, which always give correct results in metastable 
equilibrium, can give spurious results in a history dependent nonergodic ensemble. The statistical mechanics and 
thermodynamics developed here allow the derivation of expressions for multidomain ensemble averages which give the 
correct results for time nondissipative nonequilibrium ensembles. The fundamental origin of the peak in the specific 
heat near the glass transition has been unambiguously shown to be a signature of a transition from metastable 
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equilibrium to a nonergodic multi-domain ensemble. 

We have shown that the transient fluctuation relations for temperature and pressure quenches provide independent 
tests of the fundamental hypotheses used in our theory: that intra-domain populations are individually Boltzmann 
distributed, that except in the immediate vicinity of the glass transition the domains are robust with respect to small 
but finite variations in thermodynamic state variables, and that the inter-domain transition rates are negligible. 
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